Warning: The following aesthetics were dropped during statistical transformation: fill.
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
Warning: The following aesthetics were dropped during statistical transformation: fill.
ℹ This can happen when ggplot fails to infer the correct grouping structure in
the data.
ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
variable into a factor?
Show code
data.biph %>%mutate(sex =ifelse(sex ==1, "f","m")) %>%ggplot() +geom_density(aes(x = length_mm, fill =factor(sex)), alpha =0.5) +facet_wrap(~tot.age, scales ="free_y") +labs(caption ="length density by age")
mutate: converted 'sex' from integer to character (0 new NA)
Warning: Groups with fewer than two data points have been dropped.
Warning: Groups with fewer than two data points have been dropped.
Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_density()`).
Show code
# age at smoltification by spatial unitdata.biph %>%mutate(sex =ifelse(sex ==1, "f","m")) %>%filter(age.type =="both") %>%ggplot() +geom_density(aes(x = juv.age), alpha =0.3) +facet_wrap(~spat.unit) +labs(x ="age at smoltification") +xlim(-1, 5)
mutate: converted 'sex' from integer to character (0 new NA)
filter: removed 32,262 rows (41%), 47,297 rows remaining
Show code
# juvenile size at agedata.biph %>%filter(age.type =="juve.only") %>%ggplot(aes(x = juv.age, y =length_mm)) +geom_point(alpha =0.3) +geom_smooth(method = lm) +facet_wrap(~spat.unit)
filter: removed 47,297 rows (59%), 32,262 rows remaining
`geom_smooth()` using formula = 'y ~ x'
filter: removed 32,262 rows (41%), 47,297 rows remaining
mutate: converted 'sex' from integer to character (0 new NA)
Source models / Load samples
Show code
laa.samples <-readRDS(paste0(home,"/models_salmon/samples/salaa_samples_2025-12-22.RData"))$samples# # Or source model code and mcmc (load libs and data)# t <- Sys.time()# source(file = paste0(home,"/R/model devel/2-2_model_salmon_biphasic_v6.R"))# Sys.time() - t # approx 12 hours with 12 000 NUTS mcmc samples
Traces and ‘potential scale reduction factor’
Show code
node.sub <-grep("sig.l", colnames(laa.samples[[1]]), value =TRUE)laa.samples[, node.sub[], drop =FALSE] %>%mcmc_trace()
Show code
laa.samples[, node.sub[], drop =FALSE] %>%gelman.diag()
Potential scale reduction factors:
Point est. Upper C.I.
sig.l 1 1
sig.lj 1 1
Multivariate psrf
1
Show code
laa.samples[, node.sub[], drop =FALSE] %>%effectiveSize()
sig.l sig.lj
2737.001 2241.881
Show code
node.sub <-grep("^g\\[", colnames(laa.samples[[1]]), value =TRUE)laa.samples[, node.sub[], drop =FALSE] %>%mcmc_trace()
Show code
laa.samples[, node.sub[], drop =FALSE] %>%gelman.diag()
node.sub <-grep("^k.year\\[", colnames(laa.samples[[1]]), value =TRUE)laa.samples[, node.sub[sample(1:length(node.sub), 36)], drop =FALSE] %>%mcmc_trace()
Show code
pl <- laa.samples[, node.sub[], drop =FALSE] %>%gelman.diag()pl$psrf[ which(pl$psrf[, 1] >1.1),]
Point est. Upper C.I.
Show code
pl <- laa.samples[, node.sub[], drop =FALSE] %>%effectiveSize()print("effective sample size")
ungroup: no grouping variables remain
distinct: removed 79,547 rows (>99%), 12 rows remaining
Joining with `by = join_by(su)`
left_join: added one column (spat.unit)
> rows only in x 0
> rows only in data.biph %>% distinct(.. ( 0)
> matched rows 36,000
> ========
> rows total 36,000
Joining with `by = join_by(su, spat.unit)`
left_join: added one column (m.lat)
> rows only in x 0
> rows only in var.lats ( 0)
> matched rows 36,000
> ========
> rows total 36,000
mutate: new variable 'r' (factor) with 12 unique values and 0% NA
ungroup: no grouping variables remain
distinct: removed 79,547 rows (>99%), 12 rows remaining
Joining with `by = join_by(su)`
left_join: added one column (spat.unit)
> rows only in x 0
> rows only in data.biph %>% distinct(.. ( 0)
> matched rows 36,000
> ========
> rows total 36,000
Joining with `by = join_by(su, spat.unit)`
left_join: added one column (m.lat)
> rows only in x 0
> rows only in var.lats ( 0)
> matched rows 36,000
> ========
> rows total 36,000
mutate: new variable 'r' (factor) with 12 unique values and 0% NA
ungroup: no grouping variables remain
distinct: removed 79,547 rows (>99%), 12 rows remaining
Joining with `by = join_by(su)`
left_join: added one column (spat.unit)
> rows only in x 0
> rows only in data.biph %>% distinct(.. ( 0)
> matched rows 2,016,000
> ===========
> rows total 2,016,000
Joining with `by = join_by(su, spat.unit)`
left_join: added one column (m.lat)
> rows only in x 0
> rows only in var.lats ( 0)
> matched rows 2,016,000
> ===========
> rows total 2,016,000
mutate: converted 'sex' from integer to character (0 new NA)
new variable 'r' (factor) with 12 unique values and 0% NA
ungroup: no grouping variables remain
distinct: removed 79,547 rows (>99%), 12 rows remaining
Joining with `by = join_by(su)`
left_join: added one column (spat.unit)
> rows only in x 0
> rows only in data.biph %>% distinct(.. ( 0)
> matched rows 72,000
> ========
> rows total 72,000
Joining with `by = join_by(su, spat.unit)`
left_join: added one column (m.lat)
> rows only in x 0
> rows only in var.lats ( 0)
> matched rows 72,000
> ========
> rows total 72,000
mutate: converted 'sex' from integer to character (0 new NA)
new variable 'r' (factor) with 12 unique values and 0% NA
Show code
ggsave(paste0(home, "/report/salaa_linf.png"), width =17, height =14, units ="cm")
Predicted length at age curves
Show code
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 3501120 187.0 5391878 288.0 NA 5391878 288.0
Vcells 12784369 97.6 24743782 188.8 102400 20553061 156.9
Show code
# based on l.mu (age.type = "both") dfp <- laa.samples[[1]][sample(1:1500,1000), ] %>%# sample from one chain to reduce object sizespread_draws(l.mu[su,smo.age,coh,sex,gy], sep =",") %>%ungroup()
ungroup: no grouping variables remain
Show code
masu <- data.biph %>%summarise(ma =max(g.year), .by = su) %>%arrange(su)
Warning: Using one column matrices in `filter()` was deprecated in dplyr 1.1.0.
ℹ Please use one dimensional logical vectors instead.
ℹ The deprecated feature was likely used in the tidylog package.
Please report the issue at <https://github.com/elbersb/tidylog/issues>.
filter: removed 14,250,000 rows (37%), 24,750,000 rows remaining
semi_join: added no columns
> rows only in x (16,101,000)
> rows only in smcsu ( 259)
> matched rows 8,649,000
> ============
> rows total 8,649,000
filter: no rows removed
ungroup: no grouping variables remain
distinct: removed 79,547 rows (>99%), 12 rows remaining
Joining with `by = join_by(su)`
left_join: added one column (spat.unit)
> rows only in x 0
> rows only in data.biph %>% distinct(.. ( 0)
> matched rows 8,649,000
> ===========
> rows total 8,649,000
Warning: Removed 722 rows containing non-finite outside the scale range
(`stat_boxplot()`).
Show code
ggsave(paste0(home, "/report/salaa_cur.png"), width =17, height =14, units ="cm")
Warning: Removed 722 rows containing non-finite outside the scale range
(`stat_boxplot()`).
ungroup: no grouping variables remain
distinct: removed 47,285 rows (>99%), 12 rows remaining
Joining with `by = join_by(su)`
left_join: added one column (spat.unit)
> rows only in x 0
> rows only in data.b %>% distinct(su,.. ( 0)
> matched rows 2,016,000
> ===========
> rows total 2,016,000
mutate: converted 'sex' from integer to character (0 new NA)
new variable 'year' (double) with 28 unique values and 0% NA
mutate: new variable 'm_length' (double) with 56 unique values and 0% NA
filter: removed 72,000 rows (4%), 1,944,000 rows remaining
Show code
tmpK %>%ggplot() +geom_boxplot(aes(x =factor(year), y =exp(.value), color =factor(sex)),outliers =FALSE) +facet_wrap(~sex) +theme_minimal() +theme(axis.text.x =element_text(angle =90,hjust =1, size =rel(.75))) +labs(x ="year", y ="k")
Show code
tmpK %>%filter(sex =="m",!spat.unit =="Ireland_west" ) %>%ungroup() %>%ggplot() +geom_boxplot(aes(x = year, y =exp(.value), group = year), color ="#00BFC4", outliers =FALSE) +geom_line(aes(y =exp(m_length), x = year), color ="#00BFC4", alpha =0.4) +facet_wrap(~spat.unit, scales ="free") +theme_minimal() +theme(axis.text.x =element_text(angle =90,hjust =1, size =rel(.75))) +labs(x ="year", y ="k", title ="males", caption ="line is global median by year and sex")
ggsave(paste0(home, "/report/salaa_km.png"), width =17, height =14, units ="cm")tmpK %>%filter(sex =="f",!spat.unit =="Ireland_west" ) %>%ungroup() %>%ggplot() +geom_boxplot(aes(x = year, y =exp(.value), group = year), color ="#F8766D", outliers =FALSE) +geom_line(aes(y =exp(m_length), x = year), color ="#F8766D", alpha =0.5) +facet_wrap(~spat.unit, scales ="free") +theme_minimal() +theme(axis.text.x =element_text(angle =90,hjust =1, size =rel(.75))) +labs(x ="year", y ="k", title ="females", caption ="line is global median by year and sex")
ggsave(paste0(home, "/report/salaa_kf.png"), width =17, height =14, units ="cm")tmpK %>%filter(spat.unit =="Ireland_west") %>%ungroup() %>%ggplot() +geom_boxplot(aes(x = year, y =exp(.value), group = year, color = sex), outliers =FALSE) +geom_line(aes(y =exp(m_length), x = year, color = sex), alpha =0.5) +facet_wrap(~sex, scales ="free") +theme_minimal() +theme(axis.text.x =element_text(angle =90,hjust =1, size =rel(.75))) +labs(x ="year", y ="k", title ="Ireland_west", caption ="line is global median by year and sex")
mutate: new variable 'r' (integer) with 12 unique values and 0% NA
pivot_longer: reorganized (V1, V2, V3, V4, V5, …) into (c, value) [was 12x13, now 144x3]
mutate: converted 'c' from character to integer (0 new NA)
left_join: added 2 columns (spat.unit, m.lat)
> rows only in x 0
> rows only in var.lats ( 0)
> matched rows 144
> =====
> rows total 144
left_join: added 4 columns (spat.unit.x, m.lat.x, spat.unit.y, m.lat.y)
> rows only in x 0
> rows only in var.lats ( 0)
> matched rows 144
> =====
> rows total 144
mutate: converted 'r' from integer to factor (0 new NA)
converted 'c' from integer to factor (0 new NA)